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ABSTRACT 

In  this  Report  we  discuss  (1)  the  development  of  a  digital 
computer  program  to  compute  hydrofoil  loads  and  (2)  some  aspects 
of  hydrofoil  control.  The  program  computes  the  lift,  pitching 
moment,  and  flap  hinge  moment  on  a  two-dimensional  hydrofoil  with 
a  trailing  edge  flap  operating  near  a  free  surface  with  waves. 

The  computational  approach  involves  the  numerical  solution  of  an 
integral  equation  relating  an  upwash  distribution  to  a  kernel 
function  and  pressure  distribution.  The  pressure  distribution 
is  expanded  in  a  truncated  Glauert  series,  the  integration  is 
carried  out  numerically  using  a  Gaussian  quadrature,  and  the  co¬ 
efficients  of  the  Glauert  series  are  evaluated  by  a  minimum  error 
collocation  method. 

The  control  problem  investigated  involves  the  positioning 
of  a  pivoted  hydrofoil  by  means  of  a  servo-controlled  trim  tab. 

When  the  foil  is  pivoted  at  its  quarter  chord  and  control  is 
implemented  solely  by  means  of  a  servo  tab,  the  system  is  virtually 
uncontrollable.  However,  by  pivoting  the  foil  off  the  quarter 
chord  point  or  by  augmenting  the  servo  tab  with  a  servo  attached 
directly  to  the  foil,  the  system  can  be  controlled. 
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1.  INTRODUCTION 

One  of  the  principal  features  of  a  hydrofoil  craft  with 
fully  submerged  foils  is  its  ability  to  maneuver  in  a  seaway  with 
greater  isolation  from  the  sea  surface  than  virtually  every  other 
type  of  (surface)  vessel.  A  hydrofoil  craft  is  able  to  ’’platform1’ 
or  traverse  the  short  waves  of  a  rough  sea  without  undergoing 
significant  vertical  motion  and  to  ’’contour’’  over  high-amplitude 
long  wavelength  swells.  During  turning  maneuvers  the  craft  leans 
or  banks  into  the  turn  much  like  an  aircraft  or  bicycle,  thereby 
maintaining  a  shipboard  environment  in  which  the  lateral  loads 
are  minimal.  This  is  particularly  desirable  for  personnel  on 
board. 

To  control  a  hydrofoil  craft,  the  foils  must  be  continuously 
adjusted.  When  the  craft  is  travelling  a  straight  course  at  a 
constant  mean  height  in  waves,  the  orbital  motion  of  the  water 
results  in  a  time-varying  upwash  and  lift  fluctuations  on  the 
foil.  Typically,  the  foil  angle  of  attack  or  the  deflection  of 
flaps  is  varied  to  counter  the  lift  fluctuations  due  to  wave  mo¬ 
tion,  Foils  or  foil  control  surfaces  must  be  similarly  moved  to 
contour  over  waves  or  to  make  coordinated  (i.e.,  properly  banked) 
turns . 

Varying  the  foil  angle  of  attack  to  control  lift  forces  re¬ 
quires  especially  large,  heavy,  special-purpose  servomechanisms. 
These  servos  require  a  great  deal  of  power  (as  much  as  800  hp) 
and  high-pressure,  high-volume  hydraulic  lines  which  are  often 
hazardous.  In  order  primarily  to  reduce  the  servomechanism  re¬ 
quirements,  we  have  undertaken  a  study  of  the  unsteady  hydrody¬ 
namics  and  control  of  hydrofoils. 

The  four  hydrofoil  configurations  we  studied  are  illustrated  in 
Fig.  1.  Lift  on  the  basic  foil,  illustrated  as  configuration  1,  is 
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controlled  by  varying  the  angle  of  attack  of  the  entire  foil. 

The  foil  In  configuration  2  is  rigidly  attached  to  the  craft  and 
lift  is  controlled  by  deflecting  a  trailing-edge  flap.  Both  of 
these  configurations  are  presently  used  on  hydrofoil  craft.  Con¬ 
figuration  3  uses  a  trailing-edge  tab  to  rotate  the  foil  about  a 
pivot  point.  This  configuration  has  the  potential  for  very  low 
power  consumption  but  possesses  certain  inherent  controllability 
problems.  The  foil  in  configuration  4  uses  a  leading-edge  flap 
to  rotate  the  foil.  It  has  potentially  low-power  requirements 
and  ease  of  control  but  is  probably  more  susceptible  to  damage 
by  impact  with  submerged  debris . 

The  greatest  part  of  this  Report  deals  with  the  development 
of  numerical  methods  and  a  digital  computer  program  to  determine 
the  unsteady  loads  on  a  hydrofoil  with  a  flap  operating  in  waves 
in  the  vicinity  of  a  free  surface.  This  represents  an  extension 
of  an  earlier  work  [2]  which  did  not  include  the  effects  of  a 
flap.  We  do  not  explicitly  develop  results  for  a  leading-edge 
flap;  however,  this  configuration  is  simply  the  superposition  of 
a  rigid  foil  rotation  plus  a  flap  deflection.  We  also  discuss  a 
study  of  the  control  of  hydrofoils,  with  emphasis  on  the  tab- 
control  system  described  by  configuration  3. 
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2.  UNSTEADY  HYDRODYNAMICS 

2.1  Hydrodynamic  Coefficients 

Assuming  that  the  flow/foil  interaction  is  linear  (which  is 
valid  for  surface  wave  heights  that  are  small  compared  with  the 
foil  depth) ,  the  net  forces  and  moments  are  equivalent  to  the  sum 
of  the  loads  generated  by  the  upwash  on  a  foil  moving  steadily  at 
constant  speed  U  through  rough  water  plus  the  loads  caused  by  the 
heave  and  rotary  motion  of  a  foil  moving  through  calm  water.  The 
lift  and  moments  due  to  upwash  velocity  V  are  given  by 

Lv  =  CLV(k)V(k)  (la) 


MaV  =  CMaV(k)V(k) 


MgV  =  CM  v(k)V(k)  (1C) 

B 

where  most  of  the  variables  are  defined  in  Fig.  2  and  k  is  the  re¬ 
duced  frequency  (k=a)b/U),  nondimensionalized  by  the  semichord  b 
and  by  U. 

The  forces  and  moments  owing  to  the  motion  of  the  foil  are 


La(k)a+CLe(k)e+CLh(k)h 

(2a) 

!  (k)a+C  (k)e+C  (k)h 

aa  aB  ah 

(2b) 

m  (k)a+CM  (k)6+C„  (k)h. 

wgh 

(2c) 

The  purpose  of  our  hydrodynamic  analysis  is  to  evaluate  the 
twelve  functions  C.y(k)  in  Eqs .  1  and  2.  This  has  been  done  by 
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modifying  an  existing  digital  program  [i],  which  computes  the  un¬ 
steady  lift  and  moment  on  a  two-dimensional  foil  near  the  surface, 
to  include  the  effects  of  a  flap  and  the  upwash  caused  by  surface 
waves . 


2.2  Hydrodynamic  Analysis 

The  unsteady  pressure  distribution  over  the  surface  of  a  hy¬ 
drofoil  is  related  to  the  upwash  by  an  influence  function,  the 
kernel  function  [i],  which  is  a  complex  function  of  reduced  fre¬ 
quency,  depth,  and  Froude  number.  For  a  foil  which  has  its  leading 
edge  at  x  =  -1  and  its  trailing  edge  at  x  =  -»-l  as  shown  in  Fig.  3, 
this  relation  can  be  written  as : 

n  1 

V(x)  =  2?  /AP(S)K(x-S)d?  .  (3) 

Here,  AP  is  the  pressure  difference  across  the  foil,  V  is  the  up¬ 
wash  nondimensionalized  by  the  velocity  U,  and  K  is  the  kernel. 

For  a  foil  operating  at  finite  depth,  the  kernel  function 
must  satisfy  the  boundary  conditions  of  the  free  surface,  as  well 
as  those  of  the  foil  itself.  The  free  surface  is  modeled  by  in¬ 
cluding  a  virtual  image  of  the  foil  pressure  distribution  on  the 
opposite  side  of  the  surface.  During  operation  at  low  Froude  num¬ 
bers,  the  moving  pressure  disturbance  from  the  foil  causes  a  large 
enough  gravity  wave  to  affect  the  kernel  significantly.  The  com¬ 
plex  kernel  as  a  function  of  depth  and  Froude  number  for  two- 
dimensional  flow  is  developed  in  detail  by  Widnall  [ 2 ] ;  the  result 
is  given  here  in  Appendix  A.  For  a  rigid  foil  oscillating  in 
heave  or  pitch,  the  upwash  V(x)  is  known,  as  it  is  simply  the 
local  vertical  velocity  of  the  foil.  That  is: 
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heave:  V(x)  =  h  =  (h  eiwt} 

pitch:  V(x)  =  a  +a(x-a)  ^  =  (a  eiwt) [1  + ik(x-a) ]  . 

Here,  h  is  nondimensionalized  by  U  and  a  is  the  location  of  the 
foil  pivot  point  (see  Pig.  3).  With  the  kernel  and  upwash  known, 
Eq.  3  can  be  solved  for  AP(x). 

Numerical  technique  for  solving  AP 

For  attached  flow  on  the  hydrofoil,  the  Kutta  condition  at 
the  trailing  edge  must  be  satisfied.  This  boundary  condition  is 
satisfied  by  writing  the  pressure  distribution  as  a  Glauert  series 

AP  *  l  P„a  (x)  (*») 

n=l  n  n 

where 

_  cot (0/2)  n  =  1 

an  "  sin(n-l)0  9  n  2 

and 

x  =  -COS0  . 

By  combining  Eqs.  3  and  4,  the  velocity  may  be  written: 

V(x)  =  l  pn  +/  an(5)K(x-5)d£  .  (5) 

n-1  n  -1  n 

If  this  series  is  truncated  at  NOLT  terms*  and  the  integration 
carried  out,  Eq.  5  can  be  written  as 


*We  use  somewhat  cumbersome  terms  such  as  NOLT  in  this  text  to  be 
consistent  with  expressions  used  in  the  computer  program. 
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NOLT 

V(x)  =  I  P  I  U)  (6) 

n=l  n  n 


where 


+  1 

In  =  /  an(?)K(x-?)dC  • 


If,  in  the  numerical  solution  for  AP,  the  velocity  is  specified 
at  NP  collocation  points,  there  will  be  NP  equations  of  the  form 

g 

of  Eq.  6.  Using  a  notation  A{  },  where  {  }  is  a  matrix  with 
A  rows  and  B  columns,  this  set  of  equations  can  be  written 


M  1  mnolt  ”  , 

-  »  ty  o  {pi} 


(7) 


By  making  the  number  of  collocation  points,  NP,  larger  than  the 
number  of  terms  in  the  Glauert  series,  NOLT,  the  set  of  equations 
is  overdetermined  and  the  pressure  modes  can  be  solved  for  by  a 
method  which  minimizes  the  mean  square  error  of  the  resultant 
downwash  at  the  collocation  points.  The  coefficients  pn  which 
minimize  the  mean  square  error  satisfy  the  following  equation: 


a  (  NP  r  NOLT 

^LiK’-  Jip»w]  )■ 


for  n  -  1,2, •• • ,N0LT 

(8) 


The  result  [2]  is 


(P„)  ■  C(I„)T  !!„»->  Hn)T  IV} 


(9) 


where  {I}T  is  the  conjugate  transpose  of  {InK 
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The  lift  and  moment  acting  on  the  foil  can  be  found  by  in¬ 
tegrating  the  pressure  and  its  first  moment  over  the  chord  of  the 
foil.  The  lift  and  the  moment  about  the  foil  quarter-chord  in 
terms  of  the  Glauert  series  are  given  by 

CL  =  tt(Pi  +1/2  p2)  (10a) 

CM  =  l  (-P2  +P3)  •  (10b) 

Foil  with  trailing-edge  flap 

Addition  of  a  hinged  trailing-edge  flap  to  the  hydrofoil  ex¬ 
pands  the  set  of  unknown  force  and  moment  parameters  for  unsteady 
motions  to  include  the  hinge  moment  due  to  pitch  and  heave  and  the 
effect  of  flap  oscillation  on  lift,  and  moments  at  the  quarter 
chord  and  flap  hinge.  (The  geometry  of  the  flapped  hydrofoil  is 
shown  in  Fig.  2.) 

Hinge  moment 

The  moment  at  the  flap  hinge  can  be  computed  by  integrating 
the  moment  of  each  of  the  pressure  modes  over  the  flapped  portion 
of  the  foil  chord,  0c-*"ir.  The  result  is: 
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Si  =  Pi  /  cot  §  (cosO-cosO  )sinOdO 


3  0 


NOLT  tt 

+  l  p  /  sin(n-l)0(cos0  -  cosO  )sin0d0 
n=2  11  0  c 


IT 


px  /  (1  + cosO)cos0d© 
0_ 


[tt— 0  sin^O  )  sin3(0  ) 

c  c  ,  __  c 

-Tj - 15 - J  +  P2  - 3 - 


V, 


NOLT  -  Tsin(n- 

nlii  Pn  T  l~ HP 


-3)0.  sin(n+l )6  ' 

L  C 


37 


n+l 


+  \  P2  cos0c[n-ec-sinGccos0c] 


NOLT  -  [sin 

y  o-  cos©  p  — 

n=3  2  c  nL  n 


n0  sin(n-2)0 
c  c 

n-2 


]■ 


(11) 


Oscillating’ flap 

A  foil  with  an  oscillating  trailing-edge  flap  has  a  pressure 
distribution  with  a  logarithmic  singularity  at  the  flap  hinge, 
which  is  fit  by  the  pressure  series  [3]. 

J2_ 
pU* 

where 


=  i 


n=l 


+  pa 

n  n  c  c 


(12) 
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pc  “  ir  sinec 

2 

a  =  sine  ln(x-c) 
c 


B  =  amplitude  of  angular 
deflection  of  flap 


00 

and  J  p  a  is  the  Glauert  series. 
n=l  n  n 

A  relation  for  the  oscillating  flap  analogous  to  Eq.  5  is. 


oo  i  1 

v  =  I  pn  /  a  Kd?  +  p  /  a  Kd£  .  (13) 

n=l  n  -1  n  c  -1  c 


This  can  be  written: 


oo  1 

v*  =  I  P_  /  a  KdC  (14) 

n=l  n  -1  n 


with  the  ’’equivalent"  upwash  V*, 


1 

v*  =  V  -  pc  /  acKdC  ,  (15) 


which  can  be  calculated  directly  since  p  is  known.  Using  this 

V 

equivalent  upwash,  the  Glauert  pressure  modes  can  be  computed  in 
the  same  manner  as  was  used  for  the  plain  foil,  Eq.  9. 

In  order  to  find  the  equivalent  upwash,  the  integration  over 
the  chord  for  the  product  of  the  flap  hinge  pressure  mode  times 
the  kernel  must  be  performed.  This  integral, 

1  1 

/  a  Kd5  =  /  sine  lnU-c)2Kd£  ,  (16) 

-1  c  -1 
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has  a  logarithmic  singularity  at  the  hinge  point  which  can  be 
isolated  by  writing 


/  a  Kd£  =  /  (sin0K-sin0  K  )ln(£-c)  d£ 

X  ^  X  c  c* 


+  sin0  K  /  ln(C-c)  d£  . 

C  C 


Writing  the  kernel  as  K(x-£)  =  +  K(x-£),  where  K  is  a  func¬ 
tion  of  frequency,  depth,  and  Froude  number,  and  integrating  in  the 
second  term,  Eq.  17  becomes 

_/  a0KdC  =  _/  H(C)dC  +  sin0c[Kc  -  37^] 


X  -4  +  2c  In  +2  ln(l-c)2 

j.—  c 


where 


H(C)  =  (sin0[K  -  ^57]  -  sinec[: 


Kc-2¥T^  ln(5-c*).  (19) 


When  the  last  singularity  at  x=£  is  evaluated  analytically,  the 
equation  can  be  written 


1  x-e  c-e  1 

/  a  KdC  =  /  N(C)d5  +  /  N(C)d?  +  /  N(?)d? -Q(?)  for  x<c 
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1  c-e  x-e  1 

/  a  Kd?  =  /  NU)d£  +  /  N(£)d5  +  /  N(5)d£  -Q(?)  for  x  >  c 

-1  -1  c+e  x+e 

(21) 

where 

NU)  =  H(5)  +  2^y  ln(x-c)2  (22) 

Q(C)  =  -  ln(x-c)2  +  sin0c|fc  -  2ir^-c')] 

x  ^-l|  +  2c  ln(izfj  +  2  ln(l-c2)|.  (23) 

The  integrals  in  Egs .  20  and  21,  having  no  singularities,  may  be 
evaluated  by  a  simple  numerical  scheme. 

The  equivalent  upwash,  V*  as  given  in  Eq.  15,  may  now  be 
solved,  using  the  values  for  the  upwash  due  to  flap  oscillation: 

0  -1  <  x  £  c 

V  = 

-6  -B(x-c)b/U  c  <  x  <  1 

The  remainder  of  the  solution  for  the  pn  pressure  modes  is  the 
same  as  developed  in  Eqs .  4  -  9« 

The  lift  and  moment  coefficients  are  the  same  as  those  ex¬ 
pressed  in  Eqs.  10  and  11,  plus  the  effect  of  integrating  the 

flap  pressure  mode,  p  (Eq.  12),  over  the  foil.  An  equivalent 

c 

set  of  equations  is 
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CL  =  Tf(Pi  +|  p2)  +  P( 


/  a„dC 


-1 


CM  E  -  2  P( 

a 


/  ac(C+|)d? 


-1 


'M, 


-  i  l  pn  f  V  -  c  /  V? 
1  c  c 


-  i  pc[; 


ac(C  -  c)dC 


(24a) 


(24b) 


(24c) 


The  singularity  in  the  lift  coefficient  integration  can  be  iso¬ 
lated  leaving  a  nonsingular  integral  (Eq.  24a)  which  can  be  evalu¬ 
ated  numerically. 


/  a  d£  =  /sin6  ln(£-c)2d£ 
-1  c  -1 


(sin0  -  sine^ )ln(C-c)2  d^  +  sin©  /  ln(?-c)2d£ 
c  C  -1 


c  1  | 

/  +/  (sine  -  sine  ) ln(£-c)z  d£ 
-1  c+e 


+  sine„ 
c 


-4  +  2c  In 

1+c 

» 

1-0 

+  2 


ln(  1-c2  )i 


J 


(25) 


The  integral  in  Eq.  24c  for  the  flap  contribution  to  the  quarter-chord 
moment  is: 
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1  ,  1  .1 

/  ac(5+2)dC  =  /  ac(C-e)d?  +  (c  +  ±)  /  acd? 


c-e 


/  /l-?2  ln(C-c)2 (?-c)d? 

-1 


+1  , -  ,  ,  1 

+  /  /l-S2  ln(C-c)2  (C-c)d?  +  (c+±)  /  a  d?  , 

c+e  ^  -1  c 


(26) 


and  the  contribution  to  the  flap  hinge  moment  is  simply 

1  1  _ 

/  a  (C-c)d?  =  }  A-I2  ln(S-c)2  (C-c)d?  .  (27) 

c  c 

Unsteady  lift  in  waves 

Superposed  on  the  hydrofoil  steady  velocity  U  are  the  verti¬ 
cal  and  horizontal  components  of  the  orbital  motion  of  water  as¬ 
sociated  with  surface-wave  propagation.  Although  the  lift  owing 
to  the  horizontal  component  of  this  orbital  motion  has  been  shown 
by  O’Neill  to  be  present  for  foils  with  a  mean  value  of  lift 
or  camber  [4],  it  is  usually  small  and  shall  be  neglected  here. 

Let  us  evaluate  the  vertical  component  or  upwash  with  respect  to 
a  foil  moving  in  a  seaway. 

The  upwash  (nondimensionalized  by  U)  at  depth  D  associated 
with  plane  waves  of  amplitude  A,  traveling  at  an  angle  <|>  with  re¬ 
spect  to  a  stationary  xf,yf  coordinate  frame  (see  Fig.  4)  is  given 
by 
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PROPAGATING 
WAVE  TRAIN 


U 


FIG. 4  RELATION  OF  HYDROFOIL  TO  SURFACE  WAVES 
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V(x’y’t)  = 


Ao/kwS  "kw 


-k  D+i  (v'O  t  -  k  ,x’  +k  ,y') 


where  kx,  =  kweosi|> 


ky'  =  kwsln4, 

k  is  the  wavenumber,  2lr/wavelength,  and  g  is  the  gravitational 

w 

constant . 

To  account  for  ship  velocity  through  the  water,  we  may  de¬ 
fine  a  codirectional  coordinate  system  moving  at  the  speed  of  the 
ship  along  the  negative  x  axis  with  velocity  U  and  write 

<2«> 

Substituting  this  equation  into  Eq.  28  Tor  VI.M'.tl,  and  neglect, 
ing  y*  variations  gives 

A0/k7s  -kwD+i[(/V+Ukw  COS<!>)t  -  V'  COS*](30) 

V(x'  ,t)  =  - - —  e 

The  kx-cos*  term  describes  the  spatial  variation  of  upwash  which 
will  depend  on  the  angle  of  the  foil  with  respect  to  the  waves. 

The  lift  associated  with  this  upwash  is  solved  in  the  same 
manner  used  for  the  heave  and  pitch  upwash. 
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3.  COMPUTER  PROGRAMS 

Two  Fortran  II  programs  for  the  calculation  of  the  unsteady 
hydrodynamics  of  two-dimensional  hydrofoils  operating  near  a  free 
surface  are  presented.  The  first  program  computes  the  unsteady 
lift,  quarter  chord  moment,  and  flap  hinge  moment  owing  to  foil 
pitch  and  heave,  and  to  flap  oscillations.  The  second  program 
computes  the  unsteady  lift  and  the  moments  at  the  flap  hinge  and 
the  foil  quarter  chord  for  a  hydrofoil  platforming  in  waves. 

Both  of  the  programs  are  based  on  an  NSRDC  program  (see  Ref.  1 
for  a  description)  which  computes  the  unsteady  hydrodynamics  for 
a  foil  without  a  flap,  oscillating  in  pitch  and  heave.  Results 
are  computed  and  plotted  for  specific  cases  and  found  to  be  in 
good  agreement  with  theory. 

The  programs  are  written  in  Fortran  II  for  use  on  the  SDS  9^0 
research  computer  in  conjunction  with  the  Bolt  Beranek  and  Newman 
time  sharing  system.*  The  subroutines  used  are  the  same  as  those 
in  the  NSRDC  version  of  the  program  by  Widnall  except  for  changes 
necessary  in  going  from  Fortran  IV  to  II.  These  changes  include 
the  addition  of  an  arctan  function,  THETA  (X,Y)  and,  to  reduce 
programming  time  and  storage  space,  the  substitution  of  a  less 
sophisticated  version  of  the  exponential  integral  subroutine, 
EXPINT.  In  the  subroutine  for-  complex  matrix  inversion,  CMPINV, 
the  variable  field,  BLANK,  has  been  broken  down  into  its  real 
number  and  interger  components,  PIVOT  and  INDEX.  This  was  done 
because  the  SDS-9^0,  unlike  many  other  computers,  allocates  more 
space  for  the  storage  of  real  numbers  than  for  intergers.  The 
subroutine  storage  as  changed,  is  still  compatible  with  systems 


^Number  1.85,  using  the  21  January  1969  Fortran  operating  system. 
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using  the  same  amount  of  space  for  the  storage  of  real  and  inter- 
ger  numbers. 


3.1  Flap  Program 

Extending  the  NSRDC  program  to  compute  the  flap  hinge  moment, 
foil  pitch  moment,  and  lift,  owing  to  flap  deflection,  consists 
of  two  steps.  The  first  is  the  integration  of  each  term  of  the 
Glauert  series  over  the  flap  to  determine  the  moment  about  the 
hinge;  the  second  step  is  the  computation  of  the  complex  pressure 
distribution  over  the  entire  foil  due  to  flap  oscillation. 

The  integration  for  the  hinge  moment  is  done  using  the  ex¬ 
pression  derived  analytically  in  Eq.  11. 

The  computation  of  the  pressure  due  to  flap  deflection,  a 

c 

in  Eq.  12,  is  done  in  two  separate  steps.  The  first  step  is  the 
calculation  of  the  amplitude  of  the  logarithmic  pressure  term  for 
the  flap  (Eq.  12).  The  next  step  is  the  computation  of  the 
Glauert  pressure  series  due  to  the  equivalent  upwash  (Eq.  15). 

In  order  to  calculate  the  lift  and  moment  generated  by  the  flap, 
a  numerical  integration  of  the  flap  mode  must  be  done,  as  was 
seen  in  Eq.  20.  This  integration  is  done  by  a  Gaussian  quadra¬ 
ture  [5]  in  the  same  manner  as  all  the  numerical  integration  in 
the  original  NSRDC  program. 

The  program  divides  the  specified  number  of  collocation  points 
evenly  between  the  flapped  and  unflapped  portions  of  the  chord. 

In  each  portion  of  the  foil  the  points  are  distributed  sinusoidally 
in  order  to  concentrate  them  at  edges  where  the  computed  functions 
change  most  rapidly.  The  equations  governing  the  distribution  of 
collocation  points  are: 


20 


Report  No.  1970 


Bolt  Beranek  and  Newman  Inc. 


x 


i 


c-1 


cos 


TTl 

N 


for  -1  <  x  <  c 


x 


i 


c+1 

2 


cos 


tt! 

N 


for  c  <  x  <  1  . 


(31) 


3.2  Wave  Program 

To  compute  the  hydrodynamic  coefficients  of  a  hydrofoil 
operating  in  waves,  we  have  expanded  the  program  to  determine 
the  upwash  on  the  foil  as  a  function  of  the  input  parameters: 
foil  depth,  Froude  number,  heading  angle  (see  Fig.  4),  and  wave¬ 
length.  This  is  done  by  application  of  Eq,  30,  modified  to  give 
the  phase  of  the  coefficients  relative  to  the  maximum  vertical 
wave  velocity  at  a  specified  point  on  the  foil  chord  (e.g., 
quarter  chord  point).  The  computations  for  this  upwash  are 
handled  in  exactly  the  same  manner  as  those  for  heave  and  pitch 
in  the  original  program.  The  integral  for  flap  hinge  moment  as 
a  function  of  hinge  location  (Eq.  11)  has  been  added.  The  method 
of  collocation  point  location  is  the  same  as  used  for  the  flapped 
foil. 


3.3  Program  Flow  Charts  and  Print-Out  Listings 

Appendix  B  contains  functional  flow  charts  of  the  programs 
for  hydrodynamic  coefficients  of  a  foil  with  an  oscillating  flap 
and  for  a  foil  operating  in  waves.  There  has  been  an  attempt 
here  to  identify  variables  by  their  symbolic  names  in  the  program 
and  their  name"  ’om  the  theoretical  development.  Appendices  C, 
Dy  and  E  contains  listings  of  the  two  main  programs  plus  a  list¬ 
ing  of  the  subroutines. 
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3.4  Input-Output 

The  inputs  to  the  program  are  adapted  to  the  BBN  console- 
operated  time  sharing  computer  system.  Some  of  the  input  that 
was  introduced  many  times  was  stored  on  disk  files  for  ease  of 
operation.  In  a  punch  card  system,  this  input  would  probably  be 
more  easily  handled  by  reading  the  data  from  cards  for  each  run. 
The  inputs  for  the  flap  program,  which  follow  the  standard  For¬ 
tran  convention  for  interger  and  real  number  identification,  are: 

file  NAME 

(which  reads  from  a  file 

NP,  NOLT,  IM,  IFF,  RF,  FR,  D) 

The  operator  types  in  the  values  of  the  variables, 

CFLAP,  NFLP,  NLE,  NTF. 

using  FORMAT  (F  10.4,  31  10). 

The  inputs  to  the  wave  program  from  file  are  the  same  as 
those  for  the  flapped  foil.  The  typed-in  variables  are, 

WL,  THET,  FR,  PPA,  CFLAP,  NLE,  NTE 

using 

FORMAT  (5  F10.4,  214). 

The  Froude  number,  FR,  which  had  been  read  in  from  file  was  typed 
in  again  to  facilitate  easier  variation  of  this  parameter. 

Table  1  lists  and  describes  the  inputs  to  the  both  programs. 
Figures  5  and  6  are  sample  runs  of  the  FLAP  and  WAVE  programs. 

The  underlined  portions  were  typed  by  the  operator. 
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NP 

NOLT 

IM 


IFF 

RF 

FR 

D 

CFLAP 

NFLP 

NLE 

NTE 

WL 

THET 

PPA 


TABLE  1.  Input  for  Flap  and  Wave  Programs 

total  number  of  collocation  points  (even  number,  greater 
than  2  NOLT) 

number  of  terms  in  the  Glauert  Series 

determine  whether  or  not  virtual  image  on  opposite  side 
of  free  surface  is  accounted  for  in  kernel.  Set  IM  =  0 
to  exclude  free  surface  effects  of  imaged  Set  IM  =  any 
integer  other  than  zero  to  include  free  surface  effects 
of  virtual  image. 

=0:  (Pr  =  »)  Froude  effect  is  not  included  in  kernel 
t*0:  Froude  effect  is  included  in  kernel 

reduced  frequency  (k  =  [u)b]/[U0])  of  oscillation  (heave, 
pitch  or  flap) 

Froude  number,  F  =  (UQ)/(/gb) 
depth  in  semichords 

fraction  of  the  foil  occupied  by  flap 

number  of  terms  in  the  integration  for  flap  pressure 
term 

number  of  terms  in  the  integration  for  I  in  the  region 
£  <  x  (see  Eq.  6) 

number  of  terms  in  the  integration  for  In  in  the  region 
5  >  x  (see  Eq.  6) 

wavelength  in  semichords 

angle  (radians)  of  foil  heading 

coordinate  on  foil  to  which  wave  phase  is  related  (e.g., 
the  PPA  for  the  quarter  chord  is  -0.5) 


IMAGE  IFF  K  FROUDE  NUM  DEPTH  IN  SEMCHDS  NP  NOLT 
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FIG.  5  OUTPUT  FOR  SAMPLE  PROBLEM  OF  HYDROFOIL  WITH  OSCILLATING  FLAP 
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FIG.  6  OUTPUT  FOR  SAMPLE  PROBLEMS  OF  HYDROFOIL  IN  WAVES 
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3.5  Flap  Sample  Problem 

A  sample  problem  using  the  FLAP  program  for  a  hydrofoil  with 
a  25$  flap,  reduced  frequency  of  0.2,  and  no  effect  from  the  free 
surface  (IM=IFF=0)  is  shown  in  Fig.  5.  Twenty  collocation  points 
and  six  Glauert  terms  are  used,  while  eight  terms  are  used  in  the 
numerical  integration  (NLE=NTE=NFLP=8).  The  depth  and  Froude 
number  are  arbitrarily  read  in  as  1.0,  and  will  be  set  infinite 
because  we  have  specified  IM=IFF=0.  Tne  input  read  fr  ...  ne  file 
is  printed  and  identified.  The  fraction  of  the  chord  occupied  by 
the  flap  is  printed  out  after  the  data  from  the  file.  The  six 
columns  of  numbers  in  exponential  form  which  follow  are  the  com¬ 
plex  values  for  NOLT  terms  of  the  Glauert  series.  Here  row  num¬ 
ber  represents  the  order  of  the  term  and  the  columns  are  alter¬ 
nately  the  real  and  imaginary  parts  of  the  terms  for  heave,  pitch 
and  flap  oscillations,  in  that  order  from  left  to  right.  Thus, 
for  flap  oscillation  the  amplitude  of  the  tliird  term  in  the 
Glauert  series  is  P3  =  .24946  +  i. 11459.  The  lift  and  moment 
coefficients  for  the  three  cases  follow.  The  R  and  I  suffixes 
stand  for  the  real  and  imaginary  parts  respectively,  while  L 
stands  for  the  lift,  M  for  quarter  chord  moment  and  H  for  the 
hinge  moment.  The  "Comparison  with  Cornell"  is  the  heave  lift 
coefficient  divided  by  2tt  (i.e.,  the  coefficient  at  zero  fre¬ 
quency  without  free  surface  affects).  The  last  two  rows  of  out¬ 
put  are  a  summary  of  the  heave  and  then  pitch  output  in  the 
following  order: 

CLR,  CLI,  CMR,  CMI,  RD,  D,  NP,  NOLT. 
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3.6  Wave  Sample  Problem 

For  the  sample  problem  utilizing  the  WAVE  program,  we  have 
taken  a  case  with  the  foil  operating  at  a  very  high  Froude  corres¬ 
ponding  to  a  foil  speed  that  is  much  larger  than  the  wave  speed. 
This  condition  approximates  a  foil  moving  through  a  stationary 
sinusoidal  gust.  Again  the  free  surface  effects  have  been  ignored 
( IM  and  IFF  are  zero).  The  program  computes  values  in  waves  based 
upon  unity  orbital  velocity  at  the  surface.  Setting  D  =  0  makes 
the  amplitude  of  the  sinusoidal  upwash  at  the  foil  unity.  The 
wavelength  is  15.71  semi-chords  while  the  heading,  0.0,  is  directly 
into  the  seas.  The  Froude  number  is  104,  phase  is  computed  rela¬ 
tive  to  the  upwash  at  -0.5  (quarter  chord),  and  the  flap  is  25%  of 
the  chord.  The  output  format  is  the  same  as  for  the  flap  oscil¬ 
lation,  except  that  the  third  set  of  coefficients  is  for  operation 
in  waves,  and  the  wavelength  and  the  heading  are  printed  out. 


3.7  Operation 

For  a  given  computation  the  operator  must  select  the  values 
of  NOLT,  NP,  NFLP,  NILE,  and  NTE .  Running  time  is  affected  by  all 
of  these  parameters,  being  most  sensitive  to  the  first  two.  The 
minimum  number  of  Glauert  coefficients,  NOLT,  is  determined  by 
the  flap  size  and  is  critical  for  accuracy  of  the  flap  hinge 
moment.  For  flaps  which  are  shorter  than  half  of  the  chord,  this 
minimum  value  will  be 

NOLT  =  2tt/(tt-0  ). 

C 

(rounded  U£  to  an  interger). 

For  flaps  which  are  hinged  forward  of  the  mid  chord  NOLT  must  be 
at  least  three. 


27 


Report  No.  1970 


Bolt  Beranek  and  Newman  Inc. 


A  general  rule  for  NP  is  that  it  should  be  greater  than  2  *  NOLT. 

The  numerically  integrated  functions  are  fairly  smooth  and  values 
of  eight  for  NFLP,  NLE,  and  NTE  gave  good  results.  Values  up  to 
thirty-two  may  be  used. 


3.8  Results 

We  have  computed  the  twelve  hydrodynamic  functions  for  re¬ 
presentative  values  of  depth  and  Proude  number;  the  results  are 
illustrated  in  Pigs.  7  through  18.  The  functions  are  plotted 
in  the  complex  plane  with  several  values  of  reduced  frequency 
indicated  on  each  curve.  Phase  and  amplitude  may  be  found  di¬ 
rectly  from  the  plots  as  the  angle  and  magnitude  of  the  vector 
from  the  origin  to  the  point  of  interest  on  the  curve. 

Figures  8  and  9  show  pitch  and  heave  functions  for  foils 
in  waves  for  D  =  °°.  Strictly  speaking,  there  is  no  influence  of 
waves  at  infinite  depth;  this  notation  merely  means  that  the 
free  surface  effects  have  not  been  included  but  that  the  foil 
penetrates  a  sinusoidal  gust.  We  have  compared  this  solution 
with  the  Sears  function  [0]  and  have  also  compared  the  functions 
for  foil  motion  at  infinite  depth  with  Theodorsens  function  [7]. 
The  agreement  is  within  one  percent. 
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■  Fr  =  10  0=10 

ORBITAL  VELOCITY 
AT  SURFACE  *  1.0 
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•  Fr  =  oo  D  =  ® 

PHASED  TO  WAVE 

Vmox  =  1.0 

DOWNWASH  AT  1/4 

A  Fr=  10  0  =  1.0 

CHORD 

WAVE  HEIGHT  A^=1.0 

MID  CHORD  PHASE 

■  Fr=oo  D=® 

A;*  TO 
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FIG.  13  C,  PITCH  OSCILLATION 
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4.  CONTROL  PROBLEMS 

Freely  pivoted  foils,  controlled  by  trailing  edge  servo  tabs, 
exhibit  the  undesirable  feature  of  initially  developing  lift  in 
the  direction  opposite  to  that  commanded.  For  example,  when  the 
tab  is  deflected  downward,  the  lift  is  upward  until  the  foil  ro¬ 
tates  sufficiently  far  to  generate  negative  lift.  This  behavior 
leads  one  to  suspect  that  the  control  of  such  foils  may  not  be 
particularly  straightforward  and  that  there  may,  in  fact,  be  some 
fundamental  control  system  limitations  affecting  the  suitability 
of  pivoted  tab-controlled  foils  for  hydrofoil  craft.  Consequently, 
we  have  very  briefly  examined  this  problem  to  determine  its  na¬ 
ture  and  to  outline  a  possible  course  of  investigation. 

We  have  conducted  this  study  concurrently  with  the  develop¬ 
ment  of  the  hydrodynamics  program  discussed  above.  Consequently, 
the  results  of  the  program  were  not  available  for  the  control 
system  study  and  we  use  Theodorsens  equations  for  unsteady  two- 
dimensional  flow  without  free-surface  effects  [7].  Nondimension- 
alizing  the  equations  for  lift  and  foil  moment,  utilizing  Newton1 s 
laws  for  the  rigid  body  response  of  the  foil  and  ship,  and  putting 
the  equations  into  Laplace  transform  notation  gives 

-[  (m+1)  (j)+2](|)h  =  C-a4>2+2(  1-a)  0+2]a 

-it- 1  [Tj i))2  +  (TJ(-T n  ) <J>— 2T 10  ](S  (32) 

and 

-Ms+[jl  +  |  +  a2)  <J> 2 — a ( 1— 2a) <f>— 2  (a  +  |j]  a 

=  [a4»+(2a+l)]*h  +  £  { [T 7  +  ( c-a )T t ] cfj2 

+  [aTn-Tl+T8+(c-a)Tj<f>+[2aT10-T4]}|3  .  (33) 
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Here  m  is  the  ship  mass  per  unit  span  of  the  foil  divided 
by  TTpb2,  4>  is  the  Laplace  operator  nondimensionalized  by  b/U,  a 
is  the  foil  pivot  position  varying  linearly  from  -1  at  the  lead¬ 
ing  edge  to  +1  at  the  trailing  edge,  M  is  a  nondimens ional  mo- 
ment  applied  directly  to  the  foil  by  a  servomechanism,  I  is  the 
foil  moment  of  inertia  per  unit  span  divided  by  irpb4,  and  T’s 
are  functions  of  the  tab  pivot  point  defined  by  c. 

For  a  flap  comprising  1035  of  the  foil  chord,  the  coefficients 
of  3  in  Eqs.  32  and  33  may  be  approximated  by 

+  ^  U+5)  for  Eq.  (32) 

+  ^  ( <P+5 )  for  Eq.  (33) 

We  also  may  assume  that  I  is  negligible. 


4.1  Control 

The  above  equations,  in  block  diagram  form,  are  shown  in 
Fig.  19.  Lifts  and  moments  owing  to  foil  and  tab  angles  and  to 
vertical  motion  are  indicated.  In  addition,  a  feedback  gain  K  is 
shown,  representing  a  loop  in  which  a  moment  proportional  to  ver¬ 
tical  position  is  applied  directly  to  the  foil.  The  transfer 
function  for  this  system  is 


h 

(4>+5)  ^<J>2+8a<(>— 8J 

W-J\ 

|  (m+l)+ma2 

<f>4  ' 

f 

^  -  ma(l-2a)j 

<t>3 

+  [l-m(2a+l)-aK]<(.2  +  2K(l-a)*+2K }  . 


(34) 


FIG. 19  BLOCK  DIAGRAM  OF  HYDROFOIL  DYNAMICS  WITH  A  SERVO  TAB  COMPRI 
10%  OF  THE  FOIL  CHORD 
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There  are  several  interesting  aspects  of  Eq.  34.  First,  the 
quadratic  polynomial  in  the  numerator  always  has  one  positive  and 
one  negative  real  root  (i.e.,  open  loop  zeroes).  Thus  a  root 
locus  always  has  a  branch  ending  in  the  right  half  plane,  despite 
the  form  of  compensation  network  employed.  Secondly,  if  K  is 
zero,  Eq.  3 4  will  have  two  poles  at  the  origin.  This  means  thao 
an  impulse  in  3  will  result  in  a  steady  velocity  in  h.  Thus,  a 
momentary  deflection  of  the  tab  causes  a  transient  that  does  not 
end  with  the  foil  realigning  itself  with  the  inflow  vector  U. 
Instead,  the  foil  becomes  aligned  with  the  vector  sum  of  U  and 
ft,  defining  an  "effective”  angle  of  attack  of  zero.  If,  in  addi¬ 
tion  to  K  being  zero,  m  is  very  large  (corresponding  in  the  limit 
as  m  +  ®  to  a  fixed  pivot  point)  the  denominator  of  Eq.  34  be¬ 
comes 

m<Ji2  [ (g-  +  a2j  <f>2  -  a(l-2a)<|>  -  (2a+l) 

which  is  always  stable  for  a  <  -1/2  and  unstable  for  a  >  -1/2 
(corresponding  to  pivoting  the  foil  upstream  and  downstream  re¬ 
spectively  of  the  quarter  chord  point). 

If  we  assume  a  =  -/12  and  that  m  =  1  (a  good  approximation 
for  typical  hydrofoil  craft)  and  apply  Routh’s  stability  criterion 
to  the  denominator  of  Eq.  3*1,  all  roots  will  lie  in  the  left  half 
plane  (off  the  real  axis)  for  0  <  K  <  1/4.  Thus,  we  see  that  by 
choosing  an  appropriate  value  of  feedback  gain  the  open  loop 
transfer  function  can  be  made  controllable. 
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5.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  STUDY 

In  this  Report  we  have  developed  a  program  to  compute  the 
unsteady  hydrodynamic  lift  and  moment  functions  for  a  two- 
dimensional  hydrofoil  operating  near  a  free  surface,  and  have 
briefly  investigated  the  control  of  hydrofoils.  In  several  ex¬ 
ample  computations,  utilizing  the  developed  program,  we  found 
that  the  effect  of  the  free  surface  is  to  modify  the  magnitude  of 
phase  of  the  load  on  the  foil.  Lift  due  to  heave  and  pitch 
oscillations,  for  example,  decrease  with  decreasing  depth  as  one 
might  expect.  The  control  studies  indicate  certain  limitations 
in  controlling  foil  lift  by  means  of  servo  tabs. 

These  studies  are  largely  preliminary  to  the  primary  objec¬ 
tive  of  designing  hydrofoil  systems  that  have  good  response 
characteristics  and  low  power  consumption.  Based  on  these  re¬ 
sults  there  are  several  control-system  studies  that  ought  to  be 
undertaken  as  well  as  refinements  in  the  hydrodynamic  analysis. 
They  are  to: 

1.  Design  a  hydrofoil  control  system  configuration  that  maxi¬ 
mizes  bandwidth  and  minimizes  control  power.  This  involves 
the  choice  of  appropriate  feedback  loops,  compensation  func¬ 
tions,  and  foil  parameter  values. 

2.  Extend  the  unsteady  hydrodynamic  analysis  to  three  dimen¬ 
sional  hydrofoils.  This  would  be  particularly  suitable  to 
foils  of  low  aspect  ratio. 
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APPENDIX  A:  KERNEL  FUNCTION  FOR  TWO  DIMENSIONAL  FLOW 
WITH  A  FREE  SURFACE  [2] 


K  =  -  - - - e  "  I  Ei(q  )  +  iri 

2ir(x2+1ld2 )  W  1  0 


[! 


!♦¥)] 


-  IT  e’'  EK-oJ)  ♦  TiT  «  q‘  +  iti  (  1  *  J-j-L) J 

+  e  -  ^Ei(qa)  +  Hi  |l  +  J|i|]  + 


ia2  -q2 
~W 


ia 3  q3 

Tfi"  e  E1(_q3) 


ia,  q, 

+  -We  EiC-qJ 


where  a  star  denotes  complex  conjugate  and 


at  =  k[l  +  f*  +  (3  +  f_1)(4f  +  l)-35] 

a2  =  k[l  +  f"-1  -  (3  +  f‘)(4f  +  l)"3*] 

a3  =  k[l  -  f“*  +  i(3  -  f_1)(4f  -  l)-^] 

a„  =  k[l  -  f1  -  1(3  -  -  l)'3*] 

q0  =  k(2d  +  ix) 
qx  =  -Sj(2d  +  ix) 
q2  =  -s  2 ( 2d  +  ix) 
q3  =  -s 3 ( 2d  -  ix) 
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q4  =  -s4(2d  -  ix) 


f  = 


wU 

g 
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APPENDIX  B 

Flow  Charts  for  Hydrodynamic  Programs 
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FUNCTIONAL  FLOW  CHART  OF  FLAP  PROGRAM 

(using  notation  from  text  and  program) 
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(  continued) 


Singularity  Contribution 


Glauert  contribution  to 
Flap  Moment 
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APPENDIX  C 


PROGRAM  TO  COMPUTE  THE  HYDRODYNAMIC  COEFFICIENTS 
FOR  A  FLAPPED  FOIL  OPERATING  IN  WAVES 


C 


M 


.  3 

h 


C 


M 


.  » 


3h 
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C  TWO  DIMENSIONAL  PROGRAMME  FOR  FOIL  IN  WAVES 

DIMENSION  AL(  10),C1?(  10),CIM0).GH'32>,Gl<32).THETA(32).3Xf?2), 
IXOf 32).DWB(20,6)»DWl<20,6),DMR(6»6),DMI(6,6).XM00)#AM<  10)  . 
2PAM0).PB{  10),HA(  10  ),HB(  10), PR (6, 6), PI (6# 61 

3,WpR( 10),WPl( Ifl) ,WHR( 10) ,CPVC( 10) .BANGL(20),NLE1 (20) .NTE1  ( 

1*201 ,WHI( 10),WWR( 10) ,WWI( 10) ,WA( 101 ,WB(10) 

COMMON  XO,GR,GI»PO?.T,jEST,IM,IFF.D.FR.NT,S1#S2»S3fi.S3I,SUR,SUI, 

lA1.A2.A3S.A3I.AttB,  A  41 
P02=1. 57079633 
TI=- .  6 
TYPE  1?01 

1001  FOPMAT ( /%  INPUT  data  FROM  FILF:S) 

CALI  OPENR ( 5 ) 

CALL  READB ( 5, NP ) 

999  FORMAT ( 15 ) 

998  FORMAT(F10.4) 

1  CALL  READB ( 5, NOlT) 

CALL  READB ( 5. NT ) 

CALL  HEADB;5,JEST) 

CALL  READB'S.IM) 

CALL  R3ADB(5.IFF) 

CALL  READB(5,PF) 

CALL  READB(E.FR) 

call  readb(5,di 

CLOSE 

ACCEPT  996.NL.THET.FR.FPA .CFLAP.NIE.NTF 
996  FOPM AT (5Fl/. U.pltt) 

XKW=4.*P02/WL 

RF=SORT(XKW)/FP.+XKW*COS(THFT) 

WK=XKW*COS(THEt) 

VWO=EXP(-XKl<*D) 

3  FORMAT  (  515, 3F 1  '’.51 

IF(NOLT)501,501,500 
500  TYPE  2 

2  FORMAT(/62H1  IMAGE  IFF  K  FRCUDE  NUM  DEPTH  IN  SEMCHo?  NP 

1  NOLT  1 

TYPE  S.IM.IFF.RF.FR.D.NP.NOLT 
6  FORMAT(/2I5.3F10.4, 12X.2I3) 

XFLAP=1 .-2.*CFLAP 
NP2=Np/2 
PI?=4.*P02 
DO  31  J=1,NP2 
DX=C0S(J*PI2/(NP+1. )  ) 

X(.7)  =  (-(XFLAP+1.)*DX+XFLAP-1.  )/2. 

I=NP+1_J 

31  X(T)=( (-XFLAP+1 . )*DX+XFLAP+1. )/2. 

DO  32  L=1.NP 

BANGL ( L ) =ANGL  ( -X ( L ) . SQBT ( 1 ,-X f L ) *X ( L ) ) )/(2.»P02! 

NLEl (L)=NLE 
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32  NTF.  1  { I.  >  =  NTF 

TFLAP=aNGL  ( -XFL AP . SORT  (  1 . -XFI.AP*XFLAP  )  ) 

TYPE  9001.CFLAP 

9001  FORMAT (///*  FLAP  CHORD'S ,  F8 . 2 . 5X. IFR ACTION  OF  AIRFOIL  CHORD?) 
SNTf =SIH ( TFLAP ) 

TYPE  51»WL.THET 

51  FORMATf/.IWAVE  LENGTHS. F8 . 3. $  SEMI-CHORDS.  HEADING=S, F8 . 3,/ ) 
DO  100  J=1.HP 
ANRLE=BANGL( J) 

NLF=NlE1{J) 

NTF=NTS 1  { J ) 

a  FORMAT(F10. 4,2110) 

th;iax=po2»2.*angle 
X< .!)=-  COS(THMAX) 

DO  6  1=1, NLE 

THF.TA(I)=THMAX/2,*(1.-GN(I,HLE)  ) 

6  XOfI)=X(J)+COS(THETA(I)  ) 

DO  7  K=1.NOLT 
CR(K)=0.0 

7  Cl ( K ) =0 , 0 

CALL  KER2{RF,NlE) 

DO  10  1=1, NLE 

AL M ) =COS ( THETA < I) /2.) /SIN ( THETA (I)/2.) 

DO  71  N=2.  NOLT 

71  AUN)=SIN((N-1)*THETA(I)) 

CV=WN(I.NLE)*THMAX/2.*AL(2) 

DO  8  K=1,NOLT 

CR'K)=AL(K)*CW*GR(T)+CR(K) 

S  CIfK)=AL(K)*CW.SI(l)+CI(K) 

IF  (J-5)  10,9,10 

9  IF(JFST)  80,10.80 

80  TYPE  22»X0(I),GR(I1,GI(I)»CW,(AL(K),K=1.3).(CR(K).K=1,3) 

10  CONTINUE 

DO  15  1=1, NTE 

THeTA(I)=THMAX+[P02-THMAX/2.)*( l.-GN(I.NTE) )  1 

15  XO(I)=X( J)+COS(THETA(I) ) 

CALL  KER2(F.F,NTE) 

DO  25  1=1, NTE 

AL  M  )  =COS  ( '"HETA  ( I ) /2  . ) /SIN  <  THETA  <  I ) /2  .  ) 

DO  151  N=2.  NOLT 
151  AL(N)=SIN( (N-1)*THETA(I) ) 

CW=WN (I,NTF)*() 02-THMAX/2. ) *SIN ( THFTA ( I ) ) 

DO  20  K=1,"CLT 

cRf k;=al(k)*cv*gr(i)+cs(k) 

20  CI(K)=AL(K)*CW*GI(I)+CI(K)  > 

IF  (JEST)  331,25,301 
301  IF  ( J-5 )  25,21.25 

?1  TYPE  22,X0(I),GR(I),GI(I),CW,(AL(K),K=1.3),(CR(K),K=1,3) 

22  F0RMAT(//5F12. 4,/, 5512.4,/)  ' 
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25  CONTINUE 
CPVC( 1  )=TI 

DO  251  N=2.  UOlT 

251  CPVC(N)=TI*COS((N-1)*THMAX) 

DO  26  K= 1 ,  NOLI 
DWR(J»K)=CR(K)+CPVC(K) 

26  DWT(J,K)=CI(K) 

IFf J-5)100,27,1’0 

27  IF(JEST)28. 130.28 

28  TYPE  29# ( CP VC (l),I=1,3) 

29  FOR«AT(/3E20.8) 

100  CONTINUE 

DO  101  1=1. NOLT 
DO  101  K=1.NOLT 
DNS (I , K 1=0.3 

101  DMI(I,K)=0.0 

DO  110  1=1, BOLT 
DO  110  K=1 , NOLT 
DO  110  J=1 ,NP 

DMP<I,K)=DVE(J.I)*DHH(J,K)+DWI(J#I>*DWI(J,K)  +DMR(I,K1 

109  DHI(I.K)*-CWI(J»I)*DWS(J,K)+DWB(J.I)*DWI(J,K)+DMI(I,K) 

110  CONTINUE 

CALL  CMPINV(DMR,D11I,PR.PI.NOLT.INdEX1) 

IFfJEST)  111,115.111 

111  TYPE  112, ( (DHR(I,J1,DNI(I,J)#J=1, BOLT). 1=1, NOLT) 

TYPE  112,{(PR(I,J).PI(I,J).J=1.NCLT),I=1,N0LT) 

112  FORMAT(//3E20.8,/,3E20.8./) 

TYPE  112, ( (DWR( J,K) ,DWI(J,K),K=1.NOLT) ,J=1 ,NP) 

115  VPR=1.0 

VHB=1.0 

DO  150  1=1, NOLT 

WPR (I ) =0 . 0 

WPI(I)=0.0 

WHP ( I ) =0 • 0 

WHI(I)=0.0 

WWR ( I ) =0 • 0 

»WI(I)=0.0 

DO  140  J=1,NP 

VPI=(X(J)+.5)*RF 

VUpsVWD*COS ( WK* (X(J)-PPA)  ) 

VWX=-VWD*SIH(WK*(X( Jl-PPA) ) 

WHB(I)=WHR(I)+VKR.DWB(J,I) 

WHI(I)=WHI(I)-VHR*DWIf J,I) 
WPB(I)=WPR(I)+VPR*DWR(J,I)+VPI*DWI(J,I) 
WPT(I)=WPI(I)-VPR»DWI( J,I)+VPI*DWR(J,I) 

WWR(I)=WWR(I)+VWR*DWB(J,I)+VWI.DWI( J,I) 

140  WWI(I)=WWI(I)-VWR*DWI(J,I)+VWI*DWR(J,I) 

IF  (JEST)  141,150,141 

141  TYPE  1U2.WHK<I)#VHI(I).WPR(I),WPI(I) 

142  FOBNAT(/4El8.8) 
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15Ci  CONTINUE 

DO  200  1  =  1,  HOD? 

PM  I)  =0.0 
PB(1)=  1.0 
MR (I! =0.0 
H3(I}=0.0 
WRf I) =0.0 
WB tl 1=0.0 
DO  200  J= 1 , NOLT 

PA(I)=PA(I)+PPfI.J>*WPR(J)-PI(I,J;»WPl(J) 
PBf.n=pB(ii+PP(i,  ji*wpi(  ji+pki,  J:  *wPBr  j) 

HR(Il=KA(Il+P!lfI,Jl*WHR( J)-PI(I,Jl*WHI(J) 
MB'I)=HB(I)+PR(I,J)*WHI(J)+PI(I,Jl*WHR(J)  i 

w A  f 1 1 =V A ( 1 1 +  P«  ( I , J ) * WWP ( J ) -PI ( I , J 1 *WWI ( J ) 

230  WBfI)=WB(I)+PR(I,J)*HHI(.n+Pl(I,Jl*WWB(J) 

TYPE  201. (HA (II ,HB(I),pA(I),PB(I) ,WA(T},WB(I) ,1=1, NOLT) 

231  FORMAT ( /6E 1 2 . 5 1 

AL ( 1 )=P02*2.-TFLAP-SNTF 
AL(2)=P02-.5*TFLAP-,5*XFIAP*SNTF 
DO  35  N  =  3,  NOLT 

35  AI,  (N)  =  .5«(SIN(N»TFLAP)/N-SIN(N-2)*TFLAP)/(N-2) 1 
AM( 1 )=-P02+.5*TFLAP+SNTF-.5*SNTF*XFLAP 
AM(2)=(SNTF**31/3. 

AM(3)=-.5*PO2+.25*TFLAP-.0625*SINMJ.*TFLAP1 
DO  36  N=4,  NOLT 

36  AM(N)  =  .25*(SIN(N-31«TFLAP) /(N-31-riN( (N  +  1)«TFLAP)/(N+1)  ) 
DO  1102  1=1. . VOLT 

1 102  AL(I)=AM(I!-XFL#P*AMI> 

CMHR=P02/4.*(H.R(31-HA(2)  1 
CM«I=P02/U.*(HB<31-HB(2)  1 
CLHR=2.*P02*(HA( 1)+.5*HA(21) 

CLHI=2.*P02*(HB( 1)+.5*HB(2)) 

CHHP=0. 

CHt»I=0. 

DO  1201  Is-.NOLT 
CHKP.=CHHR- ,5»H A ( I )  ».AL  .  I) 

1201  CHHl=CHHI-.5*HB(I)*AL(I) 

TY®E  300 

300  FORMAT(/30HO  HEAVING  OSCILLATION  ) 

TYPE320.CLHP,CLHI,cMHB.CMHI,CHHP,CKHl 
320  F0RMAT(//$CLK=S,F8.4,$  CLI=* ,  F(i .  u ,  NX,  ICMR=J ,  F8  .4  ,  UX. 

nCMl=*,F3.U,/$cHR=i;,F8.4,UX,ICHI=».F8.tt) 

CALL  COMPRfCLHS.CLHl.RF) 

TYPE  350 

350  FOPMAT(/30H0  PITCHING  OSCILLATION  ) 

C»PR=P02/4.*(PA (3l-PA(2)  ) 

CMPI=P02/4.*(PB(3)-PB(2)  1 
CLPR=2.*P02* (PA ( 1 )+,5*PA(2) 1 
CLPI=2.*P02*(PB( 1 )+.5*PB(2) ) 

CHpR=P , 

CHPl=0 . 

DO  1202  I=1.NOLX 
CHPR=CHPR-,5*PA(I)*AL(I) 


58 


;  wvFlp „ ;  1 


WED  27-JAN-71  2  5  48PM 


PAGE  1 


1202  CHPl=CHPI~.5*PB( II* AL(I) 

TYPE  320,CLPE,CLPI.CMPR,CMPI#CHPR,CHPI 

ChwP.-P02/4  .*(Wa(3)-HA(2)  ) 

CMuI  =  po2/U.*(WM  3)-WB(2) ) 

CLWP=2.*P02*(WA( 1)+.5*WAf 2) ) 

CLWI=2.*P02* (WBf 1  )  +  .5*WB(2)) 

C  W  H  R  s  0  • 

CWHl=-?. 

HO  12^4  1  =  1#  NOLT 
CWHP=CWHR-..!3*WA(  X)*AL(I) 

1204  CWhI=CWHI-#5*V»B(I)*AL(I) 

TYoE  351 

351  F0PMAT(//SIN  WAVES. S) 

TYPE  3?0#CLWP, cLWl . CMWP , CiWI, CWHR . CWHl 
360  FOPMAT(//10H  c  «  RFAL=, 5XF 1 2 . B . 1 0H  CL  IM AG  = , 5XF 1 2 . 8 , /, 
110H  'CM  REAL®, 5XF12.8, 10H  CM  IM A0=5XF 1 2.8#/) 

TYPE  99  7#  CLH*#CLHI,cMHP.CMKI#RF#FR,D,Np.NOIT 

TYPE  947,  CLP?,CLPI#CI1PP#CMPI#RF,FR#D#NP.N0LT 

997  FOPMAT(/4F'’4,8,3F4.2#2l2) 

41  CONTINUE 
GO  TO  1 

501  CONTINUE 
STOP 
ENT) 
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APPENDIX  D:  PROGRAM  TO  COMPUTE  HYDRODYNAMIC  COEFFICIENT 
FOR  A'  FOIL  WITH  AN  OSCILLATING  FLAP. 
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C  TWO  DIMENSIONAL  programme 

DIMENSION  VFR(20),VFI(20),FKR(20).FXl(20).TT(4),WFR( 10) 

1 , WFl ( 10 ) ,FA ( 1? ) * F8 ( 1 0 ) . AM (  1 0 ) 

DIMENSION  AL( 10) ,CB( 10) »CI(10) .GR ( 32 ) »GI ( 72 ) » THETA (32) *GX ( 32) , 
1X0(32), DWB(20,fi)»DWl(20»6)*DMR(6»6)*DMI(6, 6), X(20) 
2PA(10),PB(10),HA(1S),HB(''0),PR(6»6)*PI(6»6) 

3,WPR(10),WPl(10),WHR(  10).CPVC( 10),BANGL(20),NLE,(20),NTE1! 

420) ,WHI( 10) 

COMMON  XO,GB, GI, P02 ,  T, JEST, IM.IFF, D . FR , NT. S 1 „ S 2 , S3R . S3I , SdR , SUI . 
1A1.A2.A3R,A3I.A4R,A4I 
P0?=1. 57079633 

TI=-.S 
TYPE  1001 

1001  FORMAT(/$  INFUT  data  FHOM  FILE : S ' 

CALL  OPENR ( 5 ) 

CALL  BEADB ( 5,NP ) 

999  F0RMAT(I5) 

998  FORflAT(F10.4) 

1  CALL  READB ( 5,N0LT ) 

CALL  READB(5,NT) 

CALL  READB(5,JeST) 

CALL  BEADB ( 5,IM ) 

CALL  READB ( 5, IFF ) 

CALL  READB(5,RF) 

CALL  READB(5,FR) 

CALL  READB(5,D) 

CLOSE 

3  FORMAT  ( 515, 3F 10 .5 ) 

ACCEPT  4  ,CFLAP,NFLP.NLE,NTE 
IF(NOLT)  501,501,500 
500  TYPE  2 

2  FORMAT ( /62H 1  IMAGE  IFF  K  FROUDE  NUM  DEPTH  IN  SEMCHDR  NP 

1  NOLT  ) 

TYPE  5,IM,IFF,RF,FR,D,NP,N0LT 
5  FORMAT(/2I5,3F10.4. 12X.2I3) 

XF1AP=1.-2.*CFLAP 
NP2=Np/2 
PI2=4.*P02 
'DO  31  J=1,NP2 
DX=COS( J*PI2/(NP+1 , ) ) 

X(j)=(-(XFLAP+1.)»DX+XFLAP-1.)/2. 

I=NP+1-J 

31  X ( X ) = ( (-XFLAP+1 ,)»DX+XFLAP*1. )/2. 

DO  32  L=1,NP 

BAVGL ( L )=ANGL  ( -X ( L ) , SORT( 1 . -X ( L ) .X ( L ) ) )/(2.»P02) 

NT.E1  (L)=NLE 

32  NTF 1 ( L )=NTE 

TFLAP=ANGL  (-XFLAP,SQRT( 1.-XFLAP*XFLAP)  ) 

TYPE  9001 ,CFLAP 
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0001  FORNATf///*  FLAP  CHORD*,F8.2.5X. JFBACTION  OF  AIHFOIL  CHOPDS) 
SNTF=SIN(TFLAP) 
no  15102  1*1*  HP 

1002  XO  (I ) =-COS (P02*2.*BAMGL(T) ) -XFLAP 
CALL  KEB2 ( RF»  HP ) 

DO  1003  1*1, NP 
FKI(I)=GI(I)*SNTF 

1003  FKP(I)=(GR(I)-1./(U.*P02*XO(I)) )*SNTF 

CSST=-4.*2.*ALOG( 1 .-XFLAP**2) +2. *XFLAp*ALOG( ( 1 . +XFLAP)/( 1 . -XFLAP )1 
PC*+1./(U.*P02*SMTF) 

DO  100  J=1.»P 
ANGLE=BANGL(J) 

NLF=NLE1(J) 

NTE=NTE1( J) 

4  FORMAT(F10.4,3I10) 

THMAX=P02»2. ‘ANGLE 
Xf J)=-COS(THHAX) 

DO  6  1*1, NLE 

THFTA(I)=THMAX/2.*M.-GH(I.NLE)  ) 

6  XO(I)=X(J)+COS(THETA(I) ) 

DO  7  k*1»NOLT 
CR(K)=0.0 

7  CI(K)=0.0 

CALL  KER2 ( RF, NlE ) 

DO  10  1=1, NlE 

AL(  1)=COS(THETa(I)/2.)/STN(THETA(I)/2.) 

DO  71  M=2,  NOLI 

71  Al(N)=SIN( (N-1 )*THETA(I)  ) 

CW=WH(I,NLF)*TH«AX/2.*ALf2) 

DO  8  K= 1 , NOLT 
CRfK)=Al(K)*CW*GR(I)+CR(K) 

4  CIfK)=AL(K)*CW*GI(I)+CI(K) 

IF  (J-5)  10,9,13 

9  IF ( JEST )  00,10.80 

80  TYPE  22,X0(I),GR(I),GI(I),CW,<r.L(K),K=l,3),(CR(K),K  =  1.3) 

10  continue 

DO  15  1=1, UTE 

THFTA(I)=THMAX  +  (P07-THf(AX/2, )* ' l.-GN(I.NTE) ) 

15  X0(I)=X( J)+COS(THETA(I)  ) 

CALL  KER2(RF,NTE) 

DO  25  1=1, »TE 

ALf 1)=C0S(THETA'(I)/2.)/SlH(THETA(I)/2.) 

DO  151  N=2.  NOLI 
151  AL(N)=SIN( (N-1 )*THETA(I) ) 

CW=WN(I,NTE)*(P02-THMAX/2.)*SIN(THETA(I) ) 

DO  20  K=1,N0LT 
CR(K)=AL(K)*CW*GH(I)+CR(K) 

20  CI(K)=aL(K)*CW»GI(I)+CI(K) 

IF  (JEST)  301,25,301 
301  IF  (J-5)  25,21.25 
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21  TIpE  22.XO(I),GB(I),GI(I)»CW,(AI.(K),K=1,3).(C8(Kt,K=1»3) 

22  F0pHAT(//5E12.4,/.5E12.»./) 

25  CONTINUE 
CPVC( 1 )=TX 

DO  251  N*2,  NOtT 

251  CPVC(H1=TI*C0S{ (B-1)*THMAX) 

DO  26  K=1,NOLT 
DHB(J.K)=CB(K)+CPVC(K) 

26  DHI(J.K)=CI(K) 

IF(J-5)291.27,291 

27  IF(JEST)28.291,28 

28  TYPE  29 » ( CPVC( I), Is  1,31 

29  FOFMAT(/3E20.81 

291  ALNXC=ALOG((X(J)-XFI.APn*2) 

VFB(J>=0. 

VFI(J)=0. 

TTi 11=0. 

TT(2)=THMAX 
TTf 3)=TFl*P 
TT ( 4 )=P02*2. 

IF  (TT{ 3)-TT(2) )  1011,1012.1012 

1011  VFp(J)*1. 

VFI(J1=(X(J)-XFLAP)*RF 
TT  ( 2 1 =TFLAP 

TTf  3)=THMAX 

1012  CONTINUE 

DO  1020  K=1,3 
DO  1021  I=1,NFLP 

TK»TAfI)*TT(K)+.5*(TT(K+1 )-TT(K) ) * ( 1 .-GH ( I.NFLp) ) 

1021  X0fI)=X<J)+C0S(IHETA(I1) 

CALL  KEB2(RF,NFI.P1 
DO  1022  I=1,NFLP 
SN=SIN(THETA(I) ) 

CW=WN ( I.NFLP) * .5* (TT (K+1 1 -TT(K 1 )*«N 
ALNZC=ALOG( (XFLAP+COSf THETA (1)1 1**2) 
VFR(J)=VFB(J)-PC*CW»( (S»*(GB(I1-1./(P02*<*.*X0(I) )  ) 

1-FKP.(  J)  )*ALNZC+SN/fPO2*«.*X0(in  *ALNXC) 

'022  VFI(J1=VFI(J)-PC*CW*((S»*GI(I)-FKI(J)).ALNZC) 

■020  CONTINUE 

VFP(J)=VFH(J)-PC*(..5»X(J)*ALNXC  +  FKR ( J ) »CNST ) 

VFI( J1=VFI(J)-PC*(FKI(J)»CNST) 

130  CONTINUE 

DO  101  I=1,N0LT 
DO  101  K=1.N0LT 
DNFd.K  1=0.0 
101  D.NI  ( I ,  K  )  =0 , 0 

DO  110  I=1,N0LT 
DO  110  K=1.N0LT 
DO  110  J=1,NP 

DMB(I,K)=DWP.(J.I)*DWB(J,K)+DWI(J.I)*DKI(J.K)  +D«B(I.K) 
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109  DMT(I,K)=-0«I(  J,T.)«DH'UJ.K)+DWR(J.I)*DWl(J.X)+DMI<I,K) 

life  COWTIHUE 

CALL  C«PIHV(DHR,DMI,PB.PI.HOLT.I»DEXl) 

IF(JEST)  111.115.111 

111  TYPE  112,<(DHR(I,J1,DNI(I.J),J=1,N0LT>.I=1.S0LT) 

TYPE  1 1 2 , ( (PR<I,J) .PI(I.J) ,J=1. SOLTI ,I=1.NOLT) 

112  FORHA?(//3e20.S,/,3E2C».8./) 

TYPE  112, { (DW5(J.K),DWI{J,K),K=1,N0LT1 ,J=1.HP) 

115  VPR=1.0 

VHP=1 .0 

00  150  1=1. MOLT 
MFR ( I ) =0 • 

VFI(I)=0. 

WPR(I)=0.0 

HPI(I)=0.0 

«Hfi(I)=0.0 

MHX ( X  >  =0 • 0 

00  140  J=1,NP 

VPI=(X(J)+.5)»RF 

HFR(I)=WFR(I)+VFR(.1)*DWR< J.I)+VFI(J)*DWI(J.I) 
MFI(I)=WFI(I)-VFR( J)*DWlf J.D+VFI' J)*DVR( J.I) 
WHR(I)=HHR(I)+VHR*pWR{ J.I) 

WHI(I)=WHI(I)-VHR»PWI(J,T) 

MPR(I)=WPR(I)+vPR*DME(J»X)+VPI*D«I(J»I) 

140  wpi(ii=hpi(I)-vpr*dwi(j,i)*vpi*dmr{j,i) 

IF  (JEST)  141,150,141 

141  TYPE  142.HHP(T),WHI(I).WPR(I),«PI(I) 

142  FORMAT (/4E 1 8.8) 

' 5?  costimoe 

DO  200  1=1. NOLI 
FA ( X ) =0 • 

FB (I ) =0. 

PA ( I ) =0 • 0 
PB (I ) =0 . 0 
HA ( X ) =0 . 0 
HB(I)=0.0 
PO  200  J=1.N0LT 

FB(I)=FB(I)+PP(I,J)*WFX(J)+PI(I,J1*WFR(J) 
FA(I)=FA(I)i-PR(I.J)*WF5(J)-PI(I,J)*WFI(J) 
PA(I)=PA(I)+PR(I,J)»WPR(J )-PI(I,J1*HPI(J) 
PB(I)=PB(I)+PP(I,J1*HPI(J)+PI(I,J1*WPB(J) 
HA(I)=HA(I)+PR(I, J)*WHR( J)-PI(I,J)*WHI(J) 

200  HB(I)=HB(I)+PP(I,J)*WHI( J)+PI(I.J)*HHR(J) 

TYPE  201  ,(HA(T).HB(I),PA(I).PB(I).FA(I)»FB{I).I=1.S0LT) 

201  FORMAT  ( ( /, 6E 12 • 5 ) ) 

AFL=0. 

AFH=0 . 

AFF=0. 

DO  1101  1=1 » NFLP 
THTA*.5*TFLAP»M.-GN(I.NFLP)  ) 


■;  1  •*: 
4 
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SN=SIN (THTA ) 

CS=COS(THTA) 

ALMXC=ALOG( (XFL*P+CS)**2l 
CHsWN(I,NFLP)*.5*TFLAP*SN 
AFL=AFL+CN*(SN-SNTF)*AL»XC 
AFM=AF«<+CW*SN*ALNXC*  (-CS-XFLAP1 
THTA=TFLAP+ ( P02-TFI AP* . 5 )  *  (1 .-G8 (I.NFLP ) ) 

SNsSIN(THTA) 

CSsCOS (XHTA ) 

ALNXC=ALOG( (XFLAP+CS)»*2) 

CVIsHN  (I.NFLP) *  ( P02-TFLAP*  .  5) *SN 
AFl=AFL+CS* (SN-SNTF) *AlNXC 
AFF=AFF+CW»SN*aLHXC*(-CS-XFLAP1 

1101  CONTINUE 

afl=afl+sntf*csst 

AFK=(AFH+AFF)  +  (XFL1\P+.5>*AFL 
AL( 1)=P02*2.-TFLAP-SNTF 
AL(2)=P02-.5*TFlAP-.5*XFIAP*SNTF 
DO  35  N=3,  NOLT 

35  AL(N)=.5»(SIN(n*TFLAP)/N-SIN<N-2) *TFLAP)/( N-2 )) 

AM( 1)=-P02+.5*TPLAP+SNTF-.5*SNTF*XFLAP 
AM(2)=(SNTF**3)/3. 

AM  ( 3)=-.  5*P02+.  25*'.TLAP-.0625*SIN  ( 4  .  *TFL2.P  1 
DO  36  N=4,  HOLT 

36  AM<H)  =  .2r'»fSXK{N-31»TFLAP)/(H-3)-SINUN+1)*TFLAP1/(h+1)  ) 
DO  1102  I=1,NCLX 

1102  ALfI)=A«(I)-XFLXP*AL(I) 

CMHR=P02/4.*(HA(3)-HA(2) ) 

CHHl=P02/4.*(HB(3)-HB(2) ) 

CLPR=2.*P02*(HA( 1)+.5*HA(2) ) 

CLHX=2.*P02*(HB( 1 )+.5*HB(2) ) 

CHHR=0. 

CHHI=0. 

DO  1201  1  =  1  * NOLT 
CHHR=CHHR-.5*Ha(I)*AL(I) 

1201  CHHI=CHHI-.5*HB(I)*AL(I) 

TYFE  300 

300  FOPMAT ( /30HO  HEAVING  OSCILLATION  ) 

TYPE320.CLHR. CL«I , CMHE .CMKI.CHHR.CHHI 
320  FOPMAT(//SCLH=I,F8.4,I  cli=i;,fp  .  4#  uy ,  jcmp=s  .  Fa .  u ,  ux. 

15CMl=t ,F8 ,4,/JcHR=5#  F8 .4, 4X#  $CHI=SiF8 . u ) 

CALL  COMPR(CLHB,CLHl,RF) 

TYPE  350 

350  FOPKAT(/30H?  PITCHING  OSCILLATION  ) 

CMPR=P02/4.»(PA( 3) -PM  2) ) 

CMPI=P02/4.*(PB(3)-PB(2)  ) 

CLPR=2.*P02»(PAM)  +  .5*PA(2)) 

CLPI=2.*P02*(PB( 1)+.5*PB(2)) 

CHPR=0. 

CHPl=0. 

DO  1202  I=1.NOLT 
CHPR=CHPR-.5*PA(I)*AL(I) 
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'2^2  CH?X=CHPI- . 5*PR ( I ) * AI ( T  ) 

TYPE  120,CLPF..cLPI.CK?F#CI1PI,CHPP..CHPI 
TYt>E  1301 

13/1  FOP HAT ( /  *  5  FLAP  OSCTLLATIcK i ■ 

CttFE=Pn2/U.*(FA(3)-F.M2)  >-.5*PC*AFH 
CMFl=P02/4,*(FE(3l~?F(?) \ 

CLFR  =  2.*PO?*  (Fa  (  '  )  +  .5*FA(2)  )+PC*AFL 
CLFl=2.*P02*(FR( 11+.5*FBf 2) > 

CriFH=-.S*PC*AFF 
CHFl  =  /. 

TIO  12/3  i=%nolT 
CHFP=CHFR-.S*FA(I)*AL(I) 

12/3  CHFl=CHFI-.5*FP(I)*AL(I) 

TvvE  320*  CLFP.CLFT*CMFR,C?1FI,rHFp#CHFT 

TYPE  997,  CLHB*CT  HI,CMFR.CMHI,RF,FR#P.Np#NOLT 

TYPE  QQ7,  CL©R*CLPI,CMPS*CMPI,RF*Ffl#n,HP.NOLT 

Q97  FOR«AT(/4Fl4.8.3F4.2,2r2> 

GO  TO  1 

501  CONTINUE 
STOP 
ESP 
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SUBROUTINE  CONVR(CR.CX,RF) 

PI=3. 1415193 
COR=CR/(2.*PI) 

COl*CI/(2.*PI)-RP/2. 

TYPE  10,  COB. COl 

10  FORHAK/30HO  COMPARISON  HUH  CORNELL  HR*  .P10.8.5H  HI*,F10.8) 
RETURN 
END 

SUBROUTINE  KER21RF.NCP) 

DIMENSION  GR  ( 32  )  >GK  32 )  *XO  ( 32 ) 

COMMON  X0,GH,GI,P02,X,JEST,IM,irr,D,FR,HT,S1,S2,S3R,S3J,S4R,S4X, 
1A1,A2,a3R,A3I,A4R,A4I 

CM3R(  ».,B,C,D,P,g)  =A*C.P-A*D*Q-C*B*Q-P*B*D 
CM3I( 1,B,C,D,P,g)  sA*C*Q+P*C*B+A*P*D-B*D*Q 
PI2*P02*4. 

3  DO  10  1=1, NCP 
XK=XO(X)*RF 

AB*( 1 .♦ABS  (XK)/XK)*P02 

DK*2.»'0*RF 

SINK=SIN  ( XK) 

COSK=COS  (XK) 

EXPK=EXP  (DK) 

CALL  EXPINKNT, 0.0,XK, EH, El) 

GR(I)=RF/PI2*CH3R  ( 0.0, 1 .0.COSK  ,-SINK  ,ER,EI+AB  ) 

GI( I )=RF/PI2*CM3I  (0.0,1.0,COSK  ,-SINK  ,ER,EI+AB  ) 

IF(IM)  4,10,4 

4  IF  (IFF)  6,b,6 

5  CALL  EXPINKNT,  DK,  XK,ER,EI) 

GR(I)*GR(I)-XOf  ,.)/(PI2*(XO(I)**2+4.*D**3))+RF/(2,*PI2)/2XPK 
1»CM3R  (0.0,1.0,COSK  ,-SINK  ,ER,EI+2.*AB  ) 

GI(I)=Gl(I)+RF/(2.*PI2)/BXPK  *CH3I  (0.0, 1 ,0,COSK, 

1-SINK  ,ER,EI+2.*AB  ) 

CALL  EXPINKNT,  -DK  ,XK,ER,EI) 

GR(I)=GR(I)«-BXPK  /(PI2*2,  )  *CM3R  (0.0,RF,COSK, 

1-SINK  , EK.EI) 

GI(I)=GI(I)+EXPK  /(Pl2*2.)*CM3I  (0.0,RF,COSK, 

1-SINK  , EK.EI) 

GO  TO  10 

6  CALL  EXPINKNT,  DK  ,XK,ER,EI) 

GR(I)=GR(I)-XO(i)/(P12*(XO(I)**2+4.*D**2))-RF/ ( 2, *PI2  J/EXPK 

1*CH3R  (0.0.1.0.COSK  ,-SINK  ,ER,EI+2.*AB) 

GI(I)=GI(I)-RF/(2.*PI2)/EXPX  *CM3I  (0.0, 1 ,0,COSK, 

1-SINK  ,EH,BI+2.*AB  ) 

CALL  EXPINT(NT,  -DK  ,XK,ER,EI) 

GR(I)=GR(I)-EXPK  /(PI2*2.)*CM3R  <0.0,RF,COSK, 

1-SINK  .ER.EI3 

GI(I)=GI(I)-EXPK  /( PI2*2, ) *CM3I  (0.0,RF,COSK, 

1-SINK  ,ER,EI) 

CALL  UAVBS(RF.I) 
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10  CONTINUE 
RETURN 
END 

SUBROUTINE  WAVES(RF,I) 

DIMENSION  XO(32),GR(32),GI(32) 

COMMON  X0,GR,GI,P02*X,JEST,IM,IFF,D,FR,MT,S1,S2»S3R,S3I,S4R,S4I, 
1A1»A2,A3R,A3I»A4H,A4I 
CM3H(A,B»C,D,P,U)=A*C*P-A*D*Q-C*B*Q-P*B*D 
CM3I(A,B,C,D,P,Q)sA*C*Q+P*C*B+A*P*D-B*D*0 
PI4=P02»8. 

DT=2.*D 

EX*XO(I) 

ABX=P02*2.*( 1.+ABS  (BX)/EX) 

2  FeRF»FR**2 

S1*-RF*( 1.+1./(2.*F)*( l.+SQRT  (4.*F+1,))) 

S2*-RF* ( 1.+1./(2.*F)*( 1 .-SORT  (4.*F+1,))) 

A1=RF*( 1.+1./F+(3.+1./F)*1./SQRT  (4.*F+1.)) 

A2*RP*(1.  +  l./P-(3.  +  1./F)*1./SQRT  (4,*F+1 • ) ) 

CALL  EXPINI(NT,-S1*DT,-S1*EX.ER,EI) 

GR(I)=GR(I)+EXP  (S1*DT)/PI4*CM3R  (0.0, Al,  COS  (S1*X0(I)), 

1SIN  (S1*X0(I!),ER,EI+ABX) 

GI(I)=Gl(I)+EXP  (S1*DT)/PI4*CM3I  (0.0, Al»  COS  (S1*X0(I)), 

1SIN  (S1*X0(I)),ER,EI+ABX) 

CALL  BXPINT(NT,-S2*DT,-S2*EX,ER,EI) 

GR(I)=GR(I)+:JXP  (S2*DT)/PI4*CM3R  (0.0, A2,  COS  (S2*X0(I)), 

1SIN  (S2*X0(1)),£R,EI+ABX) 

GI(I)=GI(I)+EXP  (S2»DT)/PT4*CM3I  (0.0, A2,  COS  (S2*X0(I)), 

1SIN  (S2*X0(I)),£R,!I+ABX) 

IF  (F-.25)  10,100,20 
100  STOP 
10  S-1 . -4 , *F 

RAD  =SQRT  (R) 

S3s  -RF*(  1 . - . 5/F* (  1.-RAD)  ) 

S4e  -RF*  (  1  .-•&/*’*  (  1.+RAD)  ) 

A3*  RF»( 1,-1./F-(3.-1./F)/RAD) 

A4=  RF*(1.-1./F+(3.-1./F)/RA0) 

Q3R=DT*S3 

03I=-EX*S3 

Q4R=DT*S4 

Q41s-EX*S4 

ABXM=(1.-ABS  (EX)/EX)*P02*2. 

CALL  EXPXNT<NT,g3R,03I,ER,EI) 

GR(I)»GR(I)+EXP  (-03R)/PI4*CM3R  (0.0, A3, COS  (-Q3I),SIN  (-031), 
1ER,EI+ABXM) 

GI(I)*GI(I)+EXP  (-Q3R)/PI4*CM3I  (0.0, A3, COS  (-031), SIN  (-Q3I), 
1ER,EI*ABXK) 

CALL  EXPINT(NT,04R,Q4I,ER,EI) 

GH(I)=GR(I)+EXP(-Q4R)/PI4*CM3R  (0.0,A4,COS  (-041), SIN  (-041), 
1ER,EI-ABX) 

GI(I)=Gl(I)+EXP  ( -Q4R ) /PI4*CM3I  (0.0,A4,COS  (-Q4I),SIN  (-041), 
1ER.EI-ABX) 

GO  TO  99 


69 


*  /SUBS6*6/  29  APRIL  1970  1729123 


PASS  i:2 


20  S3R=-RF*(1.-1./(2.*F)) 

S3I*-RF* ( 1 ,/ ( 2. *F) *SORT  (4.*F-1. )) 

S4R=S3B 
S4I=-S3I 
A3B=RF»( 1.-1./F) 

A3I*RF* ( 3.- 1/F ) *1  ,/SORT  (4.*F-1.) 

A4R*A3R 

A4l*-A3I 

03H=-S3R*DT-S3I*X0(I) 

03I=S3R*XO(I)-S3I*DT 

CALL  BXPINT(NT»-03H,-03I,BR,EI) 

XF  (F-.5)  21,22,22 

21  AB=P02*2.*(1.-ABS  (03D/Q3I) 

GO  TO  23 

22  AB=0, 

23  GB(I)*GR(I)+EXP  (Q3B)/PI4*CM3R  ( -A3I, A3R,COS  (Q3I5,SIR  (Q3I),BB, 

1  EX+AB) 

GI(I)=GI(I)+BXP  { Q3R)/PI4*CM3X  ( -A3I, A3R,C0S  (Q3I), SIR  (031), ER, 

1  EI*AB) 

Q42=-S4R*DT-S4J»X0(I) 

Q4l»S4R«X0(I)-S4X*D t 

CALL  EXPIHT ( NT»-Q4R,-Q4X,ER,EI ) 

IF  (F-.5)  24, 2d, 25 

24  AB=P02*2.«( 1.+ABS  (Q4I))/04l 
GO  TO  26 

25  AB=0. 

26  GI(I)=GI(I)+EXf  (Q4R)/PI4*CH3I  ( -A4I, A4R,C0S  (041), SIR  (041), ER, 

1  BI-AB) 

GR (I) =GR ( X ) +EXP  (Q4R)/PX4*CH3R  (>A4I,A4R,C0S  (041), SIR  (Q4X),ER, 

1  EX+AB) 

99  RETURN 

END 

C  COMPLEX  MATRIX  INVERSION 

SUBROUTINE  CMPINV(A,B,C,D,N,INDEX1) 

DIMENSION  A(6,6),B(6,6),C(6,6),D(6,6),PXVOT(25),XMDEX(50), 
1IPIVOT (25), SPACE (200) 

COMMON  SPACE, PIVOT, INDEX  ,IPIVOT 

C  SAVE  A  AND  INVERT  -A. 

5  M*N 

7  L*1 

10  DO  30  1*1, M 

20  DO  30  J*1,M 

30  D(I,J)*-A(I,J) 

40  CALL  MATINV  (D»M,DUMMX,0, DUMMY) 
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c  check  if  A  NAS  NON-SXNGULARc 

50  DO  70  1=1, H 

60  IF  (IPIVOT(I)-I)  250,  70,  250 

70  CONTINUE 

C  COMPUTE  C=INVERSE  OF  (A+BA(XN VERSE )B) . 

80  call  sqhult  (b,d,d,h) 

90  CALL  sqhult  (D,B,C,M) 

100  DO  128  1=1 *M 

110  DO  120  J=1,H 

120  C(I»J)=A(I*J)-C(I,J) 

130  CALL  MATINV  (C»H,DUHMY,0„ DUHMX ) 

C  CHECK  THAT  C  EXISTS. 

140  DO  170  1=1, H 

150  IF  (IPIVOT(I)-I)  155,  170,  155 

C  INVERSE  DOES  NOT  EXIST,  SET  SIGNAL  AND  RETURN. 

155  INDEX1=2 

160  RETURN 

170  CONTINUE 
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C  COMPUTE  D=-CBA(INVERSE). 

180  CALL  SQHULT  (C.D,D,H) 

220  GO  TO  (230,330),  L 

C  SUCCESSFUL  INVERSION 

230  INDEX1=1 

240  RETURN 


C  A  IS  SINGULAR,  SO  INTERCHANGE  A  AND  B  AND  TRY  AGAIN. 

250  DO  300  1=1, M 

,280  DO  308  J=1,H 

290  DUHMJ=A(I,J) 

A(I,J)=B(I,J) 

300  B(I,J)=DUHHY 

305  IF  (L-2)  310,  370,  370 

310  L=2 

GO  TO  10 

c  interchange  a  and  b,  c  and  d  WITH  CHANGED  signs. 


'  —Ml'— - 
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330  DO  350  1=1, M 

340  DO  350  J=1,M 

DUMMY=A(I,J) 

A(I,J)=B(I,J) 

B(I,J)=DUMMY 
DtlMMX  =  -C(I,J) 

C(I,J)=-D(1,J) 

350  D(I.J)=DUMMY 

360  GO  TO  230 

C  A  AND  B  BOTH  SINGULAR,  cannot  FIND  INVERSE,  SET  SIGNAL  AND  RETURN, 

370  INDEX  1  =  3 

360  RETURN 

END 

C  MATRIX  INVERSION  WITH  ACCOMPANYING  SOLUTION  OF  LINEAR  EQUATIONS 

SUBROUTINE  MATIN V ( A, N, B,M, DETERM ) 

DIMENSION  IPIVUT(25),A(6,6),B( 25,1), INDEX (25, 2), PIVOT (25), 

1  SPACE ( 200) 

COMMON  SPACE, PIVOT, INDEX, IPIVOT 

EQUIVALENCE  (IRON, GROW),  ( ICOLUM, JCOLUM ) ,  (AMAX,  T,  SNAP) 

C  INITIALIZATION 

10  DETERft=1,3 

15  ''O  20  J  =  1,N 

20  IPIVOT ( J ) =0 

30  DO  55«  1=1, N 

C  SEARCH  FOR  PIVOT  ELEMENT 

40  AMAX=3 • 0 

45  DO  105  J=1,N 

50  IF  (IPIVOT l J ) - 1 )  60,  105,  60 

60  DO  100  K=1.N 

70  IF  ( IPIVOT ( K ) -1 )  80,  100,  740 

A0  IF  (ABS  (AMAX)-ABS  (A  (J,X)))  85,  100,  100 

85  IROWsJ 

90  ICOLUM=K 

95  AMAX=A ( J,K ) 

100  CONTINUE 

105  10NTINUE 

110  IPIVOT UCOLUM)=IPIVOT(ICOLUM)+1 

C  INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 
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130  IF  (IROW-ICOLUM)  140*  260*  140 

140  DETERM3- DETERM 

160  DO  200  L=1,N 

160  SWAP=A(IROW,L) 

170  A(IROW,L)=A(ICOLUM,L) 

200  A ( ICOLUM, L )=S WAP 

205  IP(M)  260,  260*  210 

210  DO  250  L31*  M 

220  SWAP=B(IROW,L) 

230  B{IROW»L)=B(ICOLUM,L) 

250  B(IC0LUM,L)3SWAP 

260  INDEX ( I* 1 )=IROW 

270  INDEx(I,2)=ICOLUM 

310  PlVOT(I)=A(ICOLUM, ICOLUM) 

320  DETERM3DETERR*PIVOT<I) 

C  DIVIDE  PIVOT  RUN  BX  pivot  element 

330  A(ICOLUM,ICOUJfl)=1.0 

340  DO  350  L=1,N 

350  A ( ICOlUM, L ) =A ( LCOlUM, L ) /PIVOT ( I ) 

355  IF(M)  380,  380.  360 

360  DO  370  1=1, B 

370  B (ICOLUM, L)=B (ICOLUM, L) /PIVOT (I) 

C  REDUCE  NON-PIVOT  RONS 

380  DO  550  11=1, N 

390  iF(L1- ICOLUM)  400,  550,  400 

400  T=A ( L 1 , ICOLUM) 

420  A(L1,ICOLU»)=0.0 

430  DO  450  1=1, N 

450  A(L1,L)3A(l1,L)-A(ICOLUM,L)*T 

455  IF (M )  550,  550*  460 

460  DO  500  L=1,M 

500  8 (LI, L )=B ( L 1,L)-B( ICOLUM, L)*T 

550  CONTINUE 
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PACE  1 :6 


630  DO  715?  1  =  1, N 

610  LsN+1-I 

620  IF  <IHDEX(L, 1 )-INDEX(L,2) )  630,  710,  630 

630  JROK=INDEXlL,  1) 

640  JCOLUM=INOEX(L>2) 

650  DO  705  K  =  1 ,  N 

660  SH|AP=A(6,JKOW) 

670  A/K.JROW)=A(K,JCOLUM) 

700  A'K,JCOLUK)=SWAP 

705  CONTINUE 

710  CONTINUE 

740  RETURN 

END 

C  SOU  ARE  KAIBIX  MULTIPLICATION 

SUBROUTINE  SOftULT(A,B,C,N) 

DIMENSION  A(6,b),B(6,6),C(6,6),D(6,6) , COLUMN ( 25 ) , 

1  3P AC  E ( 230 ) 

COMBOS  SPACE, COLUMN 

5  M=N 

10  DO  bt.  J=1,rt 

20  DO  25  K=1,n 

25  COLUMN(K)=b(K,J) 

30  DO  50  1=1, fl 

35  C(I,J)=0.? 

U0  DO  50  K= 1 ,B 

50  C(I,J)=C{I,J)+A(I,K)*COLU«»(K) 

60  RETURN 

END 

FUNCTION  GN { 1, N ) 

DIMENSION  0(32) 

1  IF ( N-6 )  2,60,2 

2  IF ( N-8  >  3, ft, 3 

3  IF(N-lfl)  4,10,4 

4  IF  (N-16)  6,16,6 

6  IF  (N-32)  7,32,  / 

60  G{ 1 >=.9324695142 
G(2)=. 6612093865 
3(3)*. 238619161 
IF  (1-3)  61,61,62 

61  GN=G(I)  ) 

GO  TO  7  , 

62  M=7-I 

GN=-G ( M ) 

GO  TO  7 
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8  G ( 1 )=. 9602808565 
G(2)=. 79666647/4 
G(3)=. 5255324099 
G{4)=. 1834346425 
IP  (1-4)  9.9,89 

9  GM=G(I) 

GO  TO  7 

89  M=9-I 

GH“-G(») 

GO  TO  7 

10  G(1  >  =  ,9739(86529 
G(2)=, 8653633667 
G(3)=. 6799095693 
G(4)=, 43339539 
G(5)=. 1488743389 
IF  (1-5)  11,11.12 

11  GB=G(I) 

GO  TO  7 

12  11=1 1-1 
GN=-G(») 

GO  TO  7 

16  G( 1 )=. 989400935 
G(2)=. 944575 
G(3)=. 8656312 
G<4)=. 7554044 
G(5)=. 617862 
G(6)=. 45801678 
G(7)=. 2816035 

G ( 8 ) = ,3953 125 
IF  (1-8)  17,17,18 

17  GN=G(I) 

SO  TO  7 

18  11=  17-1 
SN=-G(B) 

GO  TO  7 

32  G( 1)=. 9972639 

G(2)=. 985611512 
G(3)=. 9647623 
G ( 4 ) =, 9305 
G(b)=. 89632116 
G(6)=, 849367614 
G ( 7  >  =  . 7944838 
G(8)=. 732182119 
G ( 9 )  =  ,66334426  t 
G  f  Iff  )  =  . 587715757 
G(  11  )  =  . 5(56899909 
G(12)=. 42135128 
G( 13)=. 331868602 
G(14)=, 239287362 
G(  15)  =  . 144471962 
G( 16)=. 048307666 
IF  (1-16)  33,33.34 
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33  GN=G ( I ) 

GO  TO  7 

34  (1=33-1 
GN=-G(«) 

GO  TO  7 

7  RETURN 
END 

FUNCTION  WN ( I , N  ) 
DIMENSION  N ( 32 ) 

1  IF( h-6 )  2,69,2 

2  IF  (N-8)  3.8,3 

3  IF  (N-10)  4,10,4 

4  IF( N- 16 )  6,16,6 

6  IF(N-32)  7,32,/ 

60  M( 1 )=. 1713244924 
N<2)=. 3607615731 
W(3)=. 4679139346 
IF  (1-3)  61,61,62 

61  kN— k ( I ) 

GO  TO  7 

62  M=7-I 

WN=H ( n ) 

GO  To  7 

8  W(1)=. 1912285363 
K(2)=. 2223810345 
U'<3)  =  .  3137066459 
N(4)=. 3626837834 

IF  (1-4)  9,9,89 

9  WN=W(I) 

GO  TO  7 

89  M=9-I 

MN=H(«) 

GO  TO  7 

10  4(1)=.0666713443 
W(2)=. 1494513492 
k ( 3 ) =. 2 193863625 
W/4)=. 2692667193 
W(5)=.29552U2247 
IF  (1-5)  11,11.12 

11  WN=W(I) 

GO  TO  7 
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12  0=11-1 

GO  TO  7 

16  H{ 1  )  =  . 027 152459 

•I  (2  >=3.062263924 


«'(3)  =  . 0951585117 
8(4)=. 1246289713 
«t5>=. 1495959888 
W(6)=. 1691565194 
8(7)=. 182663415 
W(8)=. 1894606105 
IF  (1-8)  1  /, 17# 18 

17  WH=»(I) 

GO  TO  7 

18  0*17-1 
W»=W(H) 

GO  TO  7 

32  W ( 1 ) =.00701 86 1 
W(2)=. 016274395 
H( 3 1=0.0253920  1 
U(4)=. 034273863 
H(5)=. 042835898 
H(6)=. 050998059 
W(7)=. 058684093 
W(8)=. 065822223 
«(9)=. 072345794 
W(10)=. 078193896 
W(11)=. 083311924 
W( 12)=. 087652093 
WF 13)=. 091173879 
K( 14)=.0938444 
W(15)=. 093872 

W( 16)=. 099654009 
IF  (1-16)  33,33,34 

33  WH=VI(I) 

50  TO  7 

34  M=33-I 
«N=W{?1) 

GO  TO  7 

7  RETURN 

END 

C  MODIFIED  PROGRAM  FOR  EXPONENTIAL  integral 
SUBROUTINE  BXPI«T(HI,X,Z,ER,EI) 
R=X*»2+X**2 

1  PI=3. 14159265 
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50  SUnB=X 
SUBI=Y 
TEBR=X 
TEMI=Y 
1 11=1.0 

S3  EHsEK-fl .0 

FACT=(EN-1.*)/(EN**2) 

TERBR=(TEBR*X-TEBI*Y)*FACT 

TERBI=(TEBI*X+T£BR*Y)*FACT 

TENH=TERBH 

T£MI=TEHBI 

SUMR=SUBR+TERBR 

SUBl=SUBI+TtRBI 

IF(R*FACT»*2-.E«S1)  60,60.53 

60  EB=.5772157+.5*AL0G(B)+SUHB 

EI=-PI+  ANGL(X,X)+SUBI 
NI=EK 
RETURN 
END 

FUNCTION  ANGL ( X,  Y ) 

C  ANGL  BETWEEN  0.  AND  2. 

IF  (Y)  10,20,30 

10  ANGL=ATAN(-X/Y)+«. 721388975 

RETURN 

20  IF ( X )  21.22,22 

21  ANGL=3.1U159265 
RETURN 

22  ANGL=0. 

RETURN 

30  ANGL=ATAN(-X/Y)+1. 57075633 

RETURN 
END 
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